maxflow.cpp

Go to the documentation of this file.
00001 /* maxflow.cpp */
00002 
00003 
00004 #include <stdio.h>
00005 #include "graph.h"
00006 #include "instances.inc"
00007 
00008 
00009 /*
00010         special constants for node->parent
00011 */
00012 #define TERMINAL ( (arc *) 1 )          /* to terminal */
00013 #define ORPHAN   ( (arc *) 2 )          /* orphan */
00014 
00015 
00016 #define INFINITE_D ((int)(((unsigned)-1)/2))            /* infinite distance to the terminal */
00017 
00018 /***********************************************************************/
00019 
00020 /*
00021         Functions for processing active list.
00022         i->next points to the next node in the list
00023         (or to i, if i is the last node in the list).
00024         If i->next is NULL iff i is not in the list.
00025 
00026         There are two queues. Active nodes are added
00027         to the end of the second queue and read from
00028         the front of the first queue. If the first queue
00029         is empty, it is replaced by the second queue
00030         (and the second queue becomes empty).
00031 */
00032 
00033 
00034 template <typename captype, typename tcaptype, typename flowtype> 
00035         inline void Graph<captype,tcaptype,flowtype>::set_active(node *i)
00036 {
00037         if (!i->next)
00038         {
00039                 /* it's not in the list yet */
00040                 if (queue_last[1]) queue_last[1] -> next = i;
00041                 else               queue_first[1]        = i;
00042                 queue_last[1] = i;
00043                 i -> next = i;
00044         }
00045 }
00046 
00047 /*
00048         Returns the next active node.
00049         If it is connected to the sink, it stays in the list,
00050         otherwise it is removed from the list
00051 */
00052 template <typename captype, typename tcaptype, typename flowtype> 
00053         inline typename Graph<captype,tcaptype,flowtype>::node* Graph<captype,tcaptype,flowtype>::next_active()
00054 {
00055         node *i;
00056 
00057         while ( 1 )
00058         {
00059                 if (!(i=queue_first[0]))
00060                 {
00061                         queue_first[0] = i = queue_first[1];
00062                         queue_last[0]  = queue_last[1];
00063                         queue_first[1] = NULL;
00064                         queue_last[1]  = NULL;
00065                         if (!i) return NULL;
00066                 }
00067 
00068                 /* remove it from the active list */
00069                 if (i->next == i) queue_first[0] = queue_last[0] = NULL;
00070                 else              queue_first[0] = i -> next;
00071                 i -> next = NULL;
00072 
00073                 /* a node in the list is active iff it has a parent */
00074                 if (i->parent) return i;
00075         }
00076 }
00077 
00078 /***********************************************************************/
00079 
00080 template <typename captype, typename tcaptype, typename flowtype> 
00081         inline void Graph<captype,tcaptype,flowtype>::set_orphan_front(node *i)
00082 {
00083         nodeptr *np;
00084         i -> parent = ORPHAN;
00085         np = nodeptr_block -> New();
00086         np -> ptr = i;
00087         np -> next = orphan_first;
00088         orphan_first = np;
00089 }
00090 
00091 template <typename captype, typename tcaptype, typename flowtype> 
00092         inline void Graph<captype,tcaptype,flowtype>::set_orphan_rear(node *i)
00093 {
00094         nodeptr *np;
00095         i -> parent = ORPHAN;
00096         np = nodeptr_block -> New();
00097         np -> ptr = i;
00098         if (orphan_last) orphan_last -> next = np;
00099         else             orphan_first        = np;
00100         orphan_last = np;
00101         np -> next = NULL;
00102 }
00103 
00104 /***********************************************************************/
00105 
00106 template <typename captype, typename tcaptype, typename flowtype> 
00107         inline void Graph<captype,tcaptype,flowtype>::add_to_changed_list(node *i)
00108 {
00109         if (changed_list && !i->is_in_changed_list)
00110         {
00111                 node_id* ptr = changed_list->New();
00112                 *ptr = (node_id)(i - nodes);
00113                 i->is_in_changed_list = true;
00114         }
00115 }
00116 
00117 /***********************************************************************/
00118 
00119 template <typename captype, typename tcaptype, typename flowtype> 
00120         void Graph<captype,tcaptype,flowtype>::maxflow_init()
00121 {
00122         node *i;
00123 
00124         queue_first[0] = queue_last[0] = NULL;
00125         queue_first[1] = queue_last[1] = NULL;
00126         orphan_first = NULL;
00127 
00128         TIME = 0;
00129 
00130         for (i=nodes; i<node_last; i++)
00131         {
00132                 i -> next = NULL;
00133                 i -> is_marked = 0;
00134                 i -> is_in_changed_list = 0;
00135                 i -> TS = TIME;
00136                 if (i->tr_cap > 0)
00137                 {
00138                         /* i is connected to the source */
00139                         i -> is_sink = 0;
00140                         i -> parent = TERMINAL;
00141                         set_active(i);
00142                         i -> DIST = 1;
00143                 }
00144                 else if (i->tr_cap < 0)
00145                 {
00146                         /* i is connected to the sink */
00147                         i -> is_sink = 1;
00148                         i -> parent = TERMINAL;
00149                         set_active(i);
00150                         i -> DIST = 1;
00151                 }
00152                 else
00153                 {
00154                         i -> parent = NULL;
00155                 }
00156         }
00157 }
00158 
00159 template <typename captype, typename tcaptype, typename flowtype> 
00160         void Graph<captype,tcaptype,flowtype>::maxflow_reuse_trees_init()
00161 {
00162         node* i;
00163         node* j;
00164         node* queue = queue_first[1];
00165         arc* a;
00166         nodeptr* np;
00167 
00168         queue_first[0] = queue_last[0] = NULL;
00169         queue_first[1] = queue_last[1] = NULL;
00170         orphan_first = orphan_last = NULL;
00171 
00172         TIME ++;
00173 
00174         while ((i=queue))
00175         {
00176                 queue = i->next;
00177                 if (queue == i) queue = NULL;
00178                 i->next = NULL;
00179                 i->is_marked = 0;
00180                 set_active(i);
00181 
00182                 if (i->tr_cap == 0)
00183                 {
00184                         if (i->parent) set_orphan_rear(i);
00185                         continue;
00186                 }
00187 
00188                 if (i->tr_cap > 0)
00189                 {
00190                         if (!i->parent || i->is_sink)
00191                         {
00192                                 i->is_sink = 0;
00193                                 for (a=i->first; a; a=a->next)
00194                                 {
00195                                         j = a->head;
00196                                         if (!j->is_marked)
00197                                         {
00198                                                 if (j->parent == a->sister) set_orphan_rear(j);
00199                                                 if (j->parent && j->is_sink && a->r_cap > 0) set_active(j);
00200                                         }
00201                                 }
00202                                 add_to_changed_list(i);
00203                         }
00204                 }
00205                 else
00206                 {
00207                         if (!i->parent || !i->is_sink)
00208                         {
00209                                 i->is_sink = 1;
00210                                 for (a=i->first; a; a=a->next)
00211                                 {
00212                                         j = a->head;
00213                                         if (!j->is_marked)
00214                                         {
00215                                                 if (j->parent == a->sister) set_orphan_rear(j);
00216                                                 if (j->parent && !j->is_sink && a->sister->r_cap > 0) set_active(j);
00217                                         }
00218                                 }
00219                                 add_to_changed_list(i);
00220                         }
00221                 }
00222                 i->parent = TERMINAL;
00223                 i -> TS = TIME;
00224                 i -> DIST = 1;
00225         }
00226 
00227         //test_consistency();
00228 
00229         /* adoption */
00230         while ((np=orphan_first))
00231         {
00232                 orphan_first = np -> next;
00233                 i = np -> ptr;
00234                 nodeptr_block -> Delete(np);
00235                 if (!orphan_first) orphan_last = NULL;
00236                 if (i->is_sink) process_sink_orphan(i);
00237                 else            process_source_orphan(i);
00238         }
00239         /* adoption end */
00240 
00241         //test_consistency();
00242 }
00243 
00244 template <typename captype, typename tcaptype, typename flowtype> 
00245         void Graph<captype,tcaptype,flowtype>::augment(arc *middle_arc)
00246 {
00247         node *i;
00248         arc *a;
00249         tcaptype bottleneck;
00250 
00251 
00252         /* 1. Finding bottleneck capacity */
00253         /* 1a - the source tree */
00254         bottleneck = middle_arc -> r_cap;
00255         for (i=middle_arc->sister->head; ; i=a->head)
00256         {
00257                 a = i -> parent;
00258                 if (a == TERMINAL) break;
00259                 if (bottleneck > a->sister->r_cap) bottleneck = a -> sister -> r_cap;
00260         }
00261         if (bottleneck > i->tr_cap) bottleneck = i -> tr_cap;
00262         /* 1b - the sink tree */
00263         for (i=middle_arc->head; ; i=a->head)
00264         {
00265                 a = i -> parent;
00266                 if (a == TERMINAL) break;
00267                 if (bottleneck > a->r_cap) bottleneck = a -> r_cap;
00268         }
00269         if (bottleneck > - i->tr_cap) bottleneck = - i -> tr_cap;
00270 
00271 
00272         /* 2. Augmenting */
00273         /* 2a - the source tree */
00274         save(middle_arc -> sister);
00275         middle_arc -> sister -> r_cap += bottleneck;
00276         save(middle_arc);
00277         middle_arc -> r_cap -= bottleneck;
00278         for (i=middle_arc->sister->head; ; i=a->head)
00279         {
00280                 a = i -> parent;
00281                 if (a == TERMINAL) break;
00282                 save(a);
00283                 a -> r_cap += bottleneck;
00284                 save(a -> sister);
00285                 a -> sister -> r_cap -= bottleneck;
00286                 if (!a->sister->r_cap)
00287                 {
00288                         set_orphan_front(i); // add i to the beginning of the adoption list
00289                 }
00290         }
00291         save(i);
00292         i -> tr_cap -= bottleneck;
00293         if (!i->tr_cap)
00294         {
00295                 set_orphan_front(i); // add i to the beginning of the adoption list
00296         }
00297         /* 2b - the sink tree */
00298         for (i=middle_arc->head; ; i=a->head)
00299         {
00300                 a = i -> parent;
00301                 if (a == TERMINAL) break;
00302                 save(a->sister);
00303                 a -> sister -> r_cap += bottleneck;
00304                 save(a);
00305                 a -> r_cap -= bottleneck;
00306                 if (!a->r_cap)
00307                 {
00308                         set_orphan_front(i); // add i to the beginning of the adoption list
00309                 }
00310         }
00311         save(i);
00312         i -> tr_cap += bottleneck;
00313         if (!i->tr_cap)
00314         {
00315                 set_orphan_front(i); // add i to the beginning of the adoption list
00316         }
00317 
00318         flow += bottleneck;
00319 }
00320 
00321 /***********************************************************************/
00322 
00323 template <typename captype, typename tcaptype, typename flowtype> 
00324         void Graph<captype,tcaptype,flowtype>::process_source_orphan(node *i)
00325 {
00326         node *j;
00327         arc *a0, *a0_min = NULL, *a;
00328         int d, d_min = INFINITE_D;
00329 
00330         /* trying to find a new parent */
00331         for (a0=i->first; a0; a0=a0->next)
00332         if (a0->sister->r_cap)
00333         {
00334                 j = a0 -> head;
00335                 if (!j->is_sink && (a=j->parent))
00336                 {
00337                         /* checking the origin of j */
00338                         d = 0;
00339                         while ( 1 )
00340                         {
00341                                 if (j->TS == TIME)
00342                                 {
00343                                         d += j -> DIST;
00344                                         break;
00345                                 }
00346                                 a = j -> parent;
00347                                 d ++;
00348                                 if (a==TERMINAL)
00349                                 {
00350                                         j -> TS = TIME;
00351                                         j -> DIST = 1;
00352                                         break;
00353                                 }
00354                                 if (a==ORPHAN) { d = INFINITE_D; break; }
00355                                 j = a -> head;
00356                         }
00357                         if (d<INFINITE_D) /* j originates from the source - done */
00358                         {
00359                                 if (d<d_min)
00360                                 {
00361                                         a0_min = a0;
00362                                         d_min = d;
00363                                 }
00364                                 /* set marks along the path */
00365                                 for (j=a0->head; j->TS!=TIME; j=j->parent->head)
00366                                 {
00367                                         j -> TS = TIME;
00368                                         j -> DIST = d --;
00369                                 }
00370                         }
00371                 }
00372         }
00373 
00374         if (i->parent = a0_min)
00375         {
00376                 i -> TS = TIME;
00377                 i -> DIST = d_min + 1;
00378         }
00379         else
00380         {
00381                 /* no parent is found */
00382                 add_to_changed_list(i);
00383 
00384                 /* process neighbors */
00385                 for (a0=i->first; a0; a0=a0->next)
00386                 {
00387                         j = a0 -> head;
00388                         if (!j->is_sink && (a=j->parent))
00389                         {
00390                                 if (a0->sister->r_cap) set_active(j);
00391                                 if (a!=TERMINAL && a!=ORPHAN && a->head==i)
00392                                 {
00393                                         set_orphan_rear(j); // add j to the end of the adoption list
00394                                 }
00395                         }
00396                 }
00397         }
00398 }
00399 
00400 template <typename captype, typename tcaptype, typename flowtype> 
00401         void Graph<captype,tcaptype,flowtype>::process_sink_orphan(node *i)
00402 {
00403         node *j;
00404         arc *a0, *a0_min = NULL, *a;
00405         int d, d_min = INFINITE_D;
00406 
00407         /* trying to find a new parent */
00408         for (a0=i->first; a0; a0=a0->next)
00409         if (a0->r_cap)
00410         {
00411                 j = a0 -> head;
00412                 if (j->is_sink && (a=j->parent))
00413                 {
00414                         /* checking the origin of j */
00415                         d = 0;
00416                         while ( 1 )
00417                         {
00418                                 if (j->TS == TIME)
00419                                 {
00420                                         d += j -> DIST;
00421                                         break;
00422                                 }
00423                                 a = j -> parent;
00424                                 d ++;
00425                                 if (a==TERMINAL)
00426                                 {
00427                                         j -> TS = TIME;
00428                                         j -> DIST = 1;
00429                                         break;
00430                                 }
00431                                 if (a==ORPHAN) { d = INFINITE_D; break; }
00432                                 j = a -> head;
00433                         }
00434                         if (d<INFINITE_D) /* j originates from the sink - done */
00435                         {
00436                                 if (d<d_min)
00437                                 {
00438                                         a0_min = a0;
00439                                         d_min = d;
00440                                 }
00441                                 /* set marks along the path */
00442                                 for (j=a0->head; j->TS!=TIME; j=j->parent->head)
00443                                 {
00444                                         j -> TS = TIME;
00445                                         j -> DIST = d --;
00446                                 }
00447                         }
00448                 }
00449         }
00450 
00451         if (i->parent = a0_min)
00452         {
00453                 i -> TS = TIME;
00454                 i -> DIST = d_min + 1;
00455         }
00456         else
00457         {
00458                 /* no parent is found */
00459                 add_to_changed_list(i);
00460 
00461                 /* process neighbors */
00462                 for (a0=i->first; a0; a0=a0->next)
00463                 {
00464                         j = a0 -> head;
00465                         if (j->is_sink && (a=j->parent))
00466                         {
00467                                 if (a0->r_cap) set_active(j);
00468                                 if (a!=TERMINAL && a!=ORPHAN && a->head==i)
00469                                 {
00470                                         set_orphan_rear(j); // add j to the end of the adoption list
00471                                 }
00472                         }
00473                 }
00474         }
00475 }
00476 
00477 /***********************************************************************/
00478 
00479 template <typename captype, typename tcaptype, typename flowtype> 
00480         flowtype Graph<captype,tcaptype,flowtype>::maxflow(bool reuse_trees, Block<node_id>* _changed_list)
00481 {
00482         node *i, *j, *current_node = NULL;
00483         arc *a;
00484         nodeptr *np, *np_next;
00485 
00486         if (!nodeptr_block)
00487         {
00488                 nodeptr_block = new DBlock<nodeptr>(NODEPTR_BLOCK_SIZE, error_function);
00489         }
00490 
00491         changed_list = _changed_list;
00492         if (maxflow_iteration == 0 && reuse_trees) { if (error_function) (*error_function)("reuse_trees cannot be used in the first call to maxflow()!"); exit(1); }
00493         if (changed_list && !reuse_trees) { if (error_function) (*error_function)("changed_list cannot be used without reuse_trees!"); exit(1); }
00494 
00495         if (reuse_trees) maxflow_reuse_trees_init();
00496         else             maxflow_init();
00497 
00498         // main loop
00499         while ( 1 )
00500         {
00501                 // test_consistency(current_node);
00502 
00503                 if ((i=current_node))
00504                 {
00505                         i -> next = NULL; /* remove active flag */
00506                         if (!i->parent) i = NULL;
00507                 }
00508                 if (!i)
00509                 {
00510                         if (!(i = next_active())) break;
00511                 }
00512 
00513                 /* growth */
00514                 if (!i->is_sink)
00515                 {
00516                         /* grow source tree */
00517                         for (a=i->first; a; a=a->next)
00518                         if (a->r_cap)
00519                         {
00520                                 j = a -> head;
00521                                 if (!j->parent)
00522                                 {
00523                                         j -> is_sink = 0;
00524                                         j -> parent = a -> sister;
00525                                         j -> TS = i -> TS;
00526                                         j -> DIST = i -> DIST + 1;
00527                                         set_active(j);
00528                                         add_to_changed_list(j);
00529                                 }
00530                                 else if (j->is_sink) break;
00531                                 else if (j->TS <= i->TS &&
00532                                          j->DIST > i->DIST)
00533                                 {
00534                                         /* heuristic - trying to make the distance from j to the source shorter */
00535                                         j -> parent = a -> sister;
00536                                         j -> TS = i -> TS;
00537                                         j -> DIST = i -> DIST + 1;
00538                                 }
00539                         }
00540                 }
00541                 else
00542                 {
00543                         /* grow sink tree */
00544                         for (a=i->first; a; a=a->next)
00545                         if (a->sister->r_cap)
00546                         {
00547                                 j = a -> head;
00548                                 if (!j->parent)
00549                                 {
00550                                         j -> is_sink = 1;
00551                                         j -> parent = a -> sister;
00552                                         j -> TS = i -> TS;
00553                                         j -> DIST = i -> DIST + 1;
00554                                         set_active(j);
00555                                         add_to_changed_list(j);
00556                                 }
00557                                 else if (!j->is_sink) { a = a -> sister; break; }
00558                                 else if (j->TS <= i->TS &&
00559                                          j->DIST > i->DIST)
00560                                 {
00561                                         /* heuristic - trying to make the distance from j to the sink shorter */
00562                                         j -> parent = a -> sister;
00563                                         j -> TS = i -> TS;
00564                                         j -> DIST = i -> DIST + 1;
00565                                 }
00566                         }
00567                 }
00568 
00569                 TIME ++;
00570 
00571                 if (a)
00572                 {
00573                         i -> next = i; /* set active flag */
00574                         current_node = i;
00575 
00576                         /* augmentation */
00577                         augment(a);
00578                         /* augmentation end */
00579 
00580                         /* adoption */
00581                         while ((np=orphan_first))
00582                         {
00583                                 np_next = np -> next;
00584                                 np -> next = NULL;
00585 
00586                                 while ((np=orphan_first))
00587                                 {
00588                                         orphan_first = np -> next;
00589                                         i = np -> ptr;
00590                                         nodeptr_block -> Delete(np);
00591                                         if (!orphan_first) orphan_last = NULL;
00592                                         if (i->is_sink) process_sink_orphan(i);
00593                                         else            process_source_orphan(i);
00594                                 }
00595 
00596                                 orphan_first = np_next;
00597                         }
00598                         /* adoption end */
00599                 }
00600                 else current_node = NULL;
00601         }
00602         // test_consistency();
00603 
00604         if (!reuse_trees || (maxflow_iteration % 64) == 0)
00605         {
00606                 delete nodeptr_block; 
00607                 nodeptr_block = NULL; 
00608         }
00609 
00610         maxflow_iteration ++;
00611         return flow;
00612 }
00613 
00614 /***********************************************************************/
00615 
00616 
00617 template <typename captype, typename tcaptype, typename flowtype> 
00618         void Graph<captype,tcaptype,flowtype>::test_consistency(node* current_node)
00619 {
00620         node *i;
00621         arc *a;
00622         int r;
00623         int num1 = 0, num2 = 0;
00624 
00625         // test whether all nodes i with i->next!=NULL are indeed in the queue
00626         for (i=nodes; i<node_last; i++)
00627         {
00628                 if (i->next || i==current_node) num1 ++;
00629         }
00630         for (r=0; r<3; r++)
00631         {
00632                 i = (r == 2) ? current_node : queue_first[r];
00633                 if (i)
00634                 for ( ; ; i=i->next)
00635                 {
00636                         num2 ++;
00637                         if (i->next == i)
00638                         {
00639                                 if (r<2){ assert(i == queue_last[r]);
00640                                 }else     assert(i == current_node);
00641                                 break;
00642                         }
00643                 }
00644         }
00645         assert(num1 == num2);
00646 
00647         for (i=nodes; i<node_last; i++)
00648         {
00649                 // test whether all edges in seach trees are non-saturated
00650                 if (i->parent == NULL) {}
00651                 else if (i->parent == ORPHAN) {}
00652                 else if (i->parent == TERMINAL)
00653                 {
00654                         if (!i->is_sink){ assert(i->tr_cap > 0);
00655                         }else             assert(i->tr_cap < 0);
00656                 }
00657                 else
00658                 {
00659                         if (!i->is_sink){ assert (i->parent->sister->r_cap > 0);
00660                         }else             assert (i->parent->r_cap > 0);
00661                 }
00662                 // test whether passive nodes in search trees have neighbors in
00663                 // a different tree through non-saturated edges
00664                 if (i->parent && !i->next)
00665                 {
00666                         if (!i->is_sink)
00667                         {
00668                                 assert(i->tr_cap >= 0);
00669                                 for (a=i->first; a; a=a->next)
00670                                 {
00671                                         if (a->r_cap > 0) assert(a->head->parent && !a->head->is_sink);
00672                                 }
00673                         }
00674                         else
00675                         {
00676                                 assert(i->tr_cap <= 0);
00677                                 for (a=i->first; a; a=a->next)
00678                                 {
00679                                         if (a->sister->r_cap > 0) assert(a->head->parent && a->head->is_sink);
00680                                 }
00681                         }
00682                 }
00683                 // test marking invariants
00684                 if (i->parent && i->parent!=ORPHAN && i->parent!=TERMINAL)
00685                 {
00686                         assert(i->TS <= i->parent->head->TS);
00687                         if (i->TS == i->parent->head->TS) assert(i->DIST > i->parent->head->DIST);
00688                 }
00689         }
00690 }

Generated on Sun Oct 26 18:22:20 2008 for maxflow-v3.0 by  doxygen 1.4.7